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ABSTRACT 



Context. Acceleration in coronal mass ejection driven shocks is currently considered the primary source of large solar energetic 
particle events. 

Aims. The solar wind, which feeds shock-accelerated particles, includes numerous ion populations, which offer much insight into 
acceleration processes. We present first simulations of shock-accelerated minor ions, in order to explore trapping dynamics and 
acceleration timescales in detail. 

Methods. We have simulated diffusive shock acceleration of minor ions ( 3 He 2+ , 4 He 2+ , 16 O s+ and 56 Fe 14+ ) and protons using a Monte 
Carlo method, where self-generated Alfvenic turbulence allows for repeated shock crossings and acceleration to high energies. 
Results. We present the effect of minor ions on wave generation, especially at low wavenumbers, and show that it is significant. We 
find that maximum ion energy is determined by the competing effects of particle escape due to focusing in an expanding flux tube 
and trapping due to the amplified turbulence. We show the dependence of cut-off energy on the particle charge to mass ratio to be 
approximately (Q/A)' 5 . 

Conclusions. We suggest that understanding the acceleration of minor ions at coronal shocks requires simulations which allow us to 
explore trapping dynamics and acceleration timescales in detail, including evolution of the turbulent trapping boundary. We conclude 
that steady-state models do not adequately describe the acceleration of heavy ions in coronal shocks. 

Key words. Acceleration of particles - Turbulence - Sun: coronal mass ejections 



1. Introduction 

Particle acceleration by coronal and interplanetary shocks driven 
by coronal mass ejections (CMEs) is widely accepted as the 
primary source of strong solar energetic particle (SEP) events. 
Particles scatter off plasma waves, cross the shock front re- 
peatedly and gain energy on each crossing. In large events, 
turbulence sufficient for extended trapping can be generated 
by the streaming of the accelerated particles themselves, as 
plasma waves in the upstream are amplified by scattering parti- 
cles. This diffusive shock acceleration mechanism, as presented 
by, e.g., IBelll dl978h . has recently been studied quant i tativel y 
as a model of the accele r ation o f SE Ps by. e.g.. | Led d2005), 
IVainio & Laitinen! (I20071 l2008h and iNg & ReameT ||2008). 
Other note worthy studies of shock-accelerated SEP events in- 
clude e.g. INg et all d!999l). ITvlka et alj d2005l) . ITvlka & Led 
(2006j) and lSandroos & Vainiol d2007l l2009albh 



Although the most abundant ion species in the solar wind, 
proton, is likely to dominate wave generation in shocks, mi- 
nor ions are scattered off the same turbulence and accelerated as 
well. If the turbulence is proton-generated, the maximum rigid- 
ity obtained by the protons determines a low-wavenumber cutoff 
in the spectrum of plasma waves. It has been suggested that this 
would prevent any ions from being accelerat ed beyond the sam e 
rigidity in a quasi-parallel shock wave 

This would then lead to a dependence of the maximum (non- 
relativistic) energy (per nucleon) of the ions on the charge to 
mass ratio of the form (Q/A) 2 . 



In this paper, we present the first simulations of particle ac- 
celeration in self-generated waves with minor ions included as 
particles contributing to the wave generation process. Instead of 
using a full spectrum of ion species in the solar wind, we limit 
ourselves to selected interesting populations, namely 3 He 2+ , 
4 He 2+ , 16 O e+ and 56 Fe I4+ . By simulating a propagating coronal 
shock for an extended period of time, we investigate the spec- 
tra and maximum attained energy of each particle population, as 
well as gauge the effect each population has on wave generation 
at different wavenumbers. 



2. Model 

In our simulation model, we numerically solve the particle trans- 
port equation through propagating representative particles using 
the guiding center approximation. We approximate the quasi- 
linear theory by employing a pitch-angle independent resonance 
condition 



fres — fci 



yv 



(1) 



where m sw is the solar wind speed, va is the Alfven speed, v is 
the particle speed, y is the Lorentz factor, fa = (l/27r)q , ,B/m,c 
is the ion cyclotron frequency, B is the magnetic flux density 
and qi = Qe and m-, = Am p are the charge and mass of the 
ion in question, and e and m v are the charge and mass of a pro- 
ton. We scatter representative particles isotropically off plasma 
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waves with scattering frequency 



2 r /res^X/res) 
V = 7T fa 



B 2 



(2) 



where P(f) is the wave power at frequency /. The initial wave 
power is scaled to give an ambient 100 keV proton mean free 
path of — 1 R Q at ro = 1.5 R Q , where R Q is the solar radius. 

Once particles have been swept up by the shock, they are 
propagated using a Monte Carlo method, focused along the mean 
magnetic field due to adiabatic invariance and traced in a su- 
perradially expanding flux tube with the solar wind speed in- 
ferred from mass conservation. For ad ditional details of our sim - 
ulation model, we refer t he r eader to Vainio & Laitinenl (120071) . 
Vainio & Laitinen (2008) and Battarbe e~al.ld2010l) . 



intense downstream turbulence (see 32.11 ). injects an unrealis- 
tically large proportion of thermal particles, we model the so- 
lar wind as consisting of two kinetic populations. The major- 
ity, 99% of particles, represents a thermal core and is passed 
directly downstream. A minority, 1% of particles, is consid- 
ered a partially sup rathermal halo and follows a ^-distribution 
dPrested et al.l 120081 This population, which has a continuous 
spectrum, is encountered by the shock in our simulation. The pa- 
rameter k receives values of 6 ... 2 going from 1.5 R e to 3.0 R e 
respectively. It should be noted, though, that the composition of 
the seed population, and its dependence on distance from the 
sun, is not known and these parameters are arbitrary. 
The average thermal speed 



w 



2Tk» 



AT- 1.5 



(4) 



2. 1 . Particle-shock interactions 

As particles encounter the propagating parallel step-profile 
shock, they scatter off the compressed plasma and may return up- 
stream with a momentum boost. The plasma compression ratio 
is solved from the Rankine-Hugoniot jump conditions. The ana- 
lytical return probability from a shock encounter for an isotropic 
particle population is given as 



Ll 2 



(3) 



where v w is the speed of the particle in the downstream plasma 
frame and u 2 the downstream plasma speed in the shock frame. 
However, the downstream-transferred populatio n is no longe r 
isotropic, unless v w » u 2 . Thus, following I Vainio et al.l (12000). 
we propagate and scatter particles in the downstream plasma 
frame up to a distance of 2A behind the shock, where A is the 
particle mean free path. At this distance, the particle has encoun- 
tered enough scatterings to warrant the assumption of isotropy 
and its representative weight is multiplied by the isotropic re- 
turn probability P Kl , sending it back towards the upstream from 
the distance of 2A with a randomized shock-bound pitch-angle. 
Particles propagate and experience small-angle scatterings in the 
downstream until they either return to the shock front or their cu- 
mulative return probability drops below 0.1 %. Simulation time 
is not advanced while the particles are in the downstream region. 
This corresponds to a situation where the downstream scattering 
is so intense that the mean residence time downstream is negligi- 
ble compared to the time between subsequent shock encounters. 
Using this assumption, we do not propagate the shape of turbu- 
lence into the downstream. 



2.2. Injected populations 

In our simulation, particles swept up by a shock propagating 
through the solar corona are accelerated and traced up along a 
grid extending to a distance of 300 R Q . In addition to a pro- 
ton population » p (r) based on the solar wind density model of 
ICranmer & van Ballegooiien (2005), we inject fully ionized he- 
lium ( 3 He 2+ and 4 He 2+ ) and partially ionized heavier elements 
( 56 Fe 14+ and 16 O s+ ) according to estimated solar wind abun- 
dance values. This results in H + -relative abundances for 4 He 2+ , 
16 6+ , 56 Fe I4+ and 3 He 2+ of 4.0 x 10~ 2 , 8.0 x 10~ 4 , 1.0 x 10~ 4 
and 1.6 x 10~ 5 , respectively. 

As a large portion of the solar wind consists of thermal parti- 
cles, and our strong step-like shock, with the assumed extremely 



is based on the radial tempe rature profile given in 
ICranmer & van Ballegooiie"nl (120051) . where T is the proton 
temperature and k% is the Boltzmann constant. All minor ion 
populations are initialised using the same distribution function 



/(v) = 



n(r)T(K + 1) 



v l n W K WT{K - 1/2) 



1 + 



(5) 



where F denotes the gamma function. 
2.3. Wave evolution 

The requirement for particle trapping in front of a shock is strong 
enough turbulence. This can be attained by wave amplification 
via the scattering of energetic particles upstream of the shock. 
A large shock-normal velocity results in large amounts of ki- 
netic energy deposited into accelerated particles and thus larger 
amounts of energy deposited into upstream waves. Also, as par- 
ticles reach higher velocities, they become resonant with waves 
of lower frequencies. 

In steady-state upstream solar wind, the evolution equation 
for a normalized wave power spectrum P = (V 2 /Bva)P can be 
written as 



dP dP ~ 8 I dP 

lH +v T r ^ P + Tf{ D "Tf 



(6) 



Here, P(r, f, f) is the Alfven wave power spectrum as a function 
of radial distance (r), frequency (/) and time (f), V = m sw + va 
is the group speed of the Alfven waves and T w is the wave 
growth rate. Dff = (V/r^f^^f' 2 ^ 3 is an ad-hoc diffusion coeffi- 
cient. This coefficient is chosen so that an unenhanced spectrum 
tends towards the form of Kolmogorov turbulence P oc /~ 5 ^ 3 at 
r e = 1 AU above the breakp oint frequency ft, = 1 mHz, as sug- 
gested by observations (e.g.. iHorburv 6130(119961) 1. Turbulence 
and diffusion magnitudes are normalized to result in an 100 kev 
proton having mean free paths of 1 R Q at r — 1.5 R Q and 54 R G 
at r = 1 AU. 

Our simulation coordinates are attached to the propagating 
coronal shock, allowing for good numerical accuracy in the near- 
shock region. Wave amplification is calculated over time inter- 
vals of 1.6 ms at grid cell boundaries. We simulate the effect 
of diffusion using a Crank-Nicholson method, and advect wave 
power along the moving grid with a Lax-Wendroff scheme uti- 
lizing a Van Leer flux limiter. 

For our simu lations, we ha ve extended the wave growth 
approximation of IVainiol (120031) to accommodate for different 
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particle masses and charges. Energy deposited into parallel- 
propagating Alfven waves per particle scattering can be written 
as AE W = -v^pAp, where p = ymv - yAm v v is the particle mo- 
mentum and Ap is the change in particle pitch-angle. Integrating 
all particles of a given ion population i at a set position gives the 
rate of change of the wave energy density as 



df 



• J"d 3 v VAyAmpV-^- F(r, v, t), 



(7) 



where F(r, v, t) is the particle distribution function. 

The pitch-angle diffusion coeffient, according to quasilinear 
theory, is 



n 1 2 n 2 ^\W(k r ) 
D m = 2k w ci (l -p ) — — 2 — 



(8) 



where w c ; = 2 w cpMy = QeB/Am p cy is the angular 
ion cyclotron frequency, k r is the resonant wavenumber and 
W(r, k, t)dk is the energy density of waves propagating parallel 
to the mean magnetic field with wavenumber in the range from 
k to k + dk. This can be written as 



cp 



2 Aym p n p (r)v' 



-(1 - p 2 )\k r \W(k r ) 



(9) 



which, because (Ap)/At = dD^/dp, yields the wave growth rate 



r (k)-±™ 
lwW ~ w dt 



j d 3 v v(l -p 2 )\k\6[k+ ^)^-- (10) 



n Q w cp 

2 n p (r) v A J \ vpj dp 



As in IVainiol (120031 we neglect the //-dependence of 
the resonance condition and replace 5(k + a> c i/vp) with 
(1/2)6 (\k\ - (L> c i/v). Now, using partial integration in p, we get 



7T Q COra 

r w (*) = T- ! 



2 M p (r) v A 



J" d 3 VV p:\k\6 (\k\-^ F (11) 



which can be further represented as 



„ , A n v x S „(r, v r , f) 

r Wii (r, v r , = - giW C p — 

2 n P (^)v A 



(12) 



where 5,, = 2ttv J j d/z vpF(r,v r ,p,t) is the particle stream- 
ing per unit velocity in the frame of the Alfven waves evalu- 
ated at the resonant particle speed v r = u)al\k\. We monitor par- 
ticle streaming by counting whenever a particle crosses a grid 
cell boundary of the shock-attached tracking grid. To correct for 
the difference of particle streaming between the frame travelling 
with the shock and the frame of the Alfven waves, we use an 
additional weighing factor of (v w ^w)/(v s jUs), where v w , p w , v s , 
and p s are the particle speed and pitch-angle in the Alfven wave 
frame and shock-attached frame, respectively. 

3. Results 

We simulate three parallel, constant velocity shocks in the 
corona starting from 1 .5 R Q . The shock-normal velocities V s are 
1250 km s~', 1500 km s~' and 1750 km s -1 . We use zero cross- 
helicity for the turbulence in the downstream, while in the up- 
stream waves propagate away from the Sun in the plasma frame. 
Turbulence is tracked on a logarithmic grid reaching out to 300 
R Q in front of the shock. 



3.1. Wave turbulence and its generation 

The left panel of Fig. fjjshows the normalized wave power spec- 
tra (multiplied with the frequency f) in front of the shock after 
160 and 640 seconds of simulation. In the right panel we show 
the wave amplification, /?;(/, z) = J T w ;(/, r s h oc k + Z,f) dt, for 
each particle population i, integrated over the whole simulation 
time. The value z is the distance from the shock front and r s h oc k 
is the position of the shock. We display wave amplification for 
the V s = 1500 km s~' run in both the measurement cell nearest 
to the shock (z = 2 ■ 10~ 5 R G ) and further out {z = 6 • 10~ 4 R Q ). 
Although minor ions are much less abundant than protons, they 
have higher charges and are resonant with lower frequencies, al- 
lowing them to generate a significant amount of turbulence close 
to the shock. S uch dominance of minor ions has been reported 
bv iLeel (fl98l for Helium ions at low frequencies in relation 
to Earth's bow shock. In our simulations, close to the shock, 
heavier ions surpass protons in wave generation at a narrow fre- 
quency range below 7 Hz. Between 7 Hz and approximately 70 
Hz, 4 He 2+ displays wave generation equal to ~ 50 % of that of 
H + . 

For the V s = 1250 km s~' shock, 4 He 2+ -powered wave gen- 
eration is equal to H + -powered generation in the region below 
100 Hz, with 4 He 2+ dominating below 30 Hz. As the shock- 
normal velocity increases, the dominance of protons on turbu- 
lence amplification increases, with the 4 He 2+ takeover moving 
to 3 Hz. In all cases, however, protons are responsible for the 
bulk of turbulence amplification, as abundant low-energy pro- 
tons generate a great deal of turbulence above 100 Hz. 

3.2. Accelerated particle populations 

In Fig. [2] we demonstrate the evolution of two particle popula- 
tions, H + and 56 Fe I4+ . Particles within the expanding flux tube 
are efficiently accelerated up to a maximum energy at which the 
turbulence can no longer trap the particles, and instead the focus- 
ing effect of the diverging magnetic field allows them to escape. 

The focused diffusi on model of particl e transport has been 
examined in detail bv lKocharov et alJdl996l) . Here, we start with 
Parker's equation, which in the fixed frame, upstream of our 
coronal shock reads 



dt dr 



3Wdr { ^ V) dp 



1 d ( dfo' 

^d?r D -d r - l 



(13) 



where fo is the isotropic part of the distribution function, D = 
(1 /3)Av is the spatial diffusion coefficient, A = v/vis the particle 
mean free path and J\ is the flux-tube cross-sectional area related 
to the focusing length L by L l = dJi/dr. Parker's equa- 
tion can be expressed, using the linear density n = d 2 N/drdp = 
Anp 2 J\fo, as a Fokker-Planck equation 



dn d 
~dt + ~dr 



( V+ l) n 



d_ 

dp 



p (dV V 
3\-d? + L r " 



d I dn) 



(14) 



which shows that the effect of focusing in the particle motion is 
two-fold: it contributes to the advection velocity by hf oc — D/L 
and to the adiabatic energy changes by pf oc = -(p/3) V/L. The 
addition to the advection velocity at large distances from the Sun 
is large, since there the waves have not yet grown to make D 
small. It is clear that particles will, on average, move away from 
the shock in the upstream region in areas where V + D/L > V s . 
This facilitates the escape of particles from the shock to the up- 
stream and the distance where V + D/L — V s can be regarded as 
the boundary of the turbulent trapping region ahead of the shock. 
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Fig. 1. Left: Evolution of wave power spectra. Right: Wave amplification factor T Wi ; integrated over 640 seconds of simulation, at 
2 • 1(T 5 R e and 6 • 1(T 4 R from the shock, where V s = 1500 km s" 1 
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Fig. 2. Particle populations for H + and 56 Fe along with the critical contour where the focusing velocity V/L exceeds the shock 
velocity. Results are shown for the simulations where V s = 1500 km s _1 , with colour contours at one magnitude intervals. A radial 
3-cell boxcar smooth function has been applied. 



This boundary, displayed in Fig. [2] as a dashed line, outlines an 
escaping population further away from the shock. In addition, 
the energy at which the boundary intersects the shock surface is 
representative of the maximum energy that the particles can be 
accelerated to at a given time. 

In the latter stages of the simulation, we find decreased wave 
generation due to lower particle densities and thus less swept-up 
particles. This, along with wave diffusion, results in the turbulent 
trapping boundary at the shock moving to lower energies, which 



causes high energy particles to escape instead of experiencing 
further acceleration. 

3.2.1. Spectral indices 

The accelerated particle populations were integrated over the 
whole upstream. We found the spectral index a for the power- 
law part of the particle spectrum by fitting a line to a chosen sec- 
tion of the log-log representation of data points. Observing the 
one magnitude contours in Fig. [2] we see that the spatial distribu- 



4 



M. Battarbee et al.: Heavy-ion Acceleration and Self-generated Waves in Coronal Shocks 



-a 

G 



0.5 



0.0 



- -0.5 



o 
(X 



-1.0 



-1.5 



Vs= 


1750 


km s -1 


S6p e 14+ 


X 


Vs= 


1500 


km s" ] 


' 56p e 14+ 


□ 


Vs= 


1250 


km s -1 


' 56p e 14+ 


A 




Vs= 


1750 


km s -1 


H + 


O 


Vs= 


1500 


km s -1 


H + 




v s = 


1250 


km s -1 


H + 


+ 



100 200 



300 
Time 



400 
["] 



500 600 



700 



Fig. 4. Temporal evolution of spectral index a for H + and 56 Fe 14+ 
populations, as fitted to the power-law section of the particle 
spectra. 

tions of iron and protons differ for both the escaping population 
and particle populations within the turbulent trapping boundary. 
For iron, this results in m uch harder p article spectra than what 
the steady-state model of iBelll £1978) suggests. Spectra along 
with more complete parameter fits are exemplified in Fig. [3] 

As the simulation continues the spectrum for heavy ions soft- 
ens. Figure|4]displays the evolution of the spectral indices a for 
H + and 56 Fe 14+ . Protons, being the dominant particle population, 
do not exhibit a softening of the spectral index, whereas the ef- 
fect is exceedingly prominent in the case of iron accelerated by 
a V s = 1250 kms -1 shock. 



3.2.2. Attained maximum energies 

When attempting to gauge the maximum energy attained by a 
particle population, we attempted to fit a power law with an ex- 
ponential cut-off to the simulated spectrum. First, we found the 
spectral index a for the power-law part of the particle spectrum, 
as in jj3.2.1l We then used this as the basis for fitting an expo- 
nential cut-off energy E c . The form used is 



where C is a fitted constant and e is chosen to fit the sharpness 
of the cut-off. In our work, we used values of e = 4 ... 2.5, with 
the value decreasing over simulation time. Figure|5]displays how 
the cut-off energy follows a ratio of mass-to-charge to the power 
of 1.5-1.6, where the exponent is significantly smaller than the 
theoretical estimation of 2. 



4. Discussion & Conclusions 

Having studied the acceleration of multiple particle populations 
through self-generated turbulence with three different coronal 
shocks, we find that during early phases of the acceleration pro- 
cess, very hard, even flat spectra can be seen for high-mass ions. 
At all but the highest frequencies, the effect of minor ions on 
wave generation is non-negligible, especially in the region di- 
rectly in front of the shock. As the shock-normal velocity in- 
creases, the deduced spectra become harder and the maximum 
energy attained increases. It is also seen that the maximum en- 
ergy dependence (Q/Af does not exhibit /3 = 2, as suggested 
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Fig. 5. Ratios of particle cut-off energy to the charge/mass num- 
ber as a log-log-plot. 



bv lZank et alj (120071) . Rather, the behaviour of cut-off energies 
is between p - 1.5 and (3 — 1 .6. 

At high energies, accelerated particles stream away from 
the shock due to focusing, as scattering particles supply insuf- 
ficient wave amplification power to trap high-energy particles to 
the shock. This causes the turbulent trapping boundary to ap- 
proach and intersect the shock at the maximum ion energy. As 
the shock-normal velocity increases, the particle spectra become 
harder and the energy at which the turbulent trapping boundary 
intersects the shock increases. 

To gauge the effect of focusing versus trapping on parti- 
cle energy, we note that out of V + D/L = V s only the diffu- 
sion coefficient D oc v 2 / (/d/res-P(/res)) depends on the parti- 
cle species. Thus, at the maximum energy edge of the trapping 
boundary, D must be same for all particles. If at low frequencies 
the wave spectrum increases from ambient to amplified levels 
with a power law P oc f b , as shown in Fig.[TJ we find 



D oc 



(!) 



-6-2 



,.6+3 



(15) 



This gives the (non-relativistic) cutoff energy per nucleon a de- 
pendency of (g/A) 2(i+2) / (i+3) . For a sharp cutoff, i.e. a purely 
rigidity-limited case, this results in a ((2/A) 2 -dependence for the 
cutoff energy, while a smoother transition results in a signifi- 
cantly weaker dependence. Due to the weak dependence of /3 on 
b, and the dynamic evolution of turbulence, care should be taken 
when inferring the turbulence spectrum shape from ion cutoff 
energies. 

To examine particle acceleration dynamics, we can calculate 
the time tr required to accelerate a particle from an injection 
rigidity Rq to a given rigidity R = p/q-, oc Av/Q, where R » /Jo- 
Assuming zero residence time in the downstream, we find 



t r =C 



f R dR' 
Jro ~K 



D(R') = Qc(G(R)-G(R )), 



(16) 



where C is a constant and G(R) is a function of rigidity based 
on the shape of the wave power spectrum. If P(f) oc /* where b 
is constant, the acceleration time to a given rigidity R is directly 
proportional to Q/A. Another item of interest is the time required 
to accelerate a heavy ion from injection speed vq to the maximum 



speed v ion = (Q/A) 



(6+2)/(6+3) 



where v p is the proton speed at 
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Fig. 3. Particle spectra after 80 seconds (left panel) and 640 seconds (right panel) of simulation, where V s = 1500 km s '. A 
power-law and an exponentional cut-off has been fitted to each population, ignoring the enhancement at the lowest energies. 



the turbulent trapping boundary. This gives 

^"trap 



(!) 



(v'f + W 



(sr- 



-b-2 



(17) 



which, using previous assumptions for b, yields an acceleration 
time independent of the charge-to-mass ratio. Thus, it is clear 
that minor ions are accelerated to the maximum rigidity of pro- 
tons much faster than the protons themselves, after which the 
ions slowly continue to gain energy until they reach the turbu- 
lent trapping boundary. This results in minor ions gaining harder, 
even fiat spectra, especially in early phases of the simulation. In 
latter phases of the simulation, the value of b increases at low 
frequencies, which leads to an increase in minor ion accelera- 
tion time and thus softer minor ion spectra. 

At later stages of the simulations, wave amplification rates 
decay in line with the decay of injection efficiency. Thus, accel- 
erated particles can stream away from the shock at lower ener- 
gies, and further acceleration to higher energies ceases. Particles 
reaching the turbulent trapping boundary propagate in space, 
forming a plateau which is not consistent with Bell's steady-state 
result. 

In conclusion, diffusive shock acceleration of protons and 
minor ions cannot realistically be represented by a steady- state 
approximation, but instead requires numerical simulations to re- 
veal the full dynamics of the acceleration process and the various 
particle populations. 
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